Dyadic Representational Similarity Analysis (DyRSA)¶
The following code implements the method described in:
Rhoads, S. A. & Marsh, A. A. (2023). Modeling the minds of friends: Dorsomedial prefrontal cortex encodes accurate dyadic meta-evaluations about how others perceive us. PsyArXiv. https://doi.org/10.31234/osf.io/g9e3z
In [1]:
%matplotlib inline
%reload_ext autoreload
%autoreload 2
In [ ]:
import numpy as np, pickle
import matplotlib.pyplot as plt, seaborn as sns
from scipy.stats import spearmanr
from urllib.request import urlopen
import nibabel as nib
from nilearn.image import new_img_like
from nilearn.plotting import view_img_on_surf
from nltools.stats import fdr
from joblib import Parallel, delayed
from tqdm import tqdm
from copy import deepcopy
import dyrsa
In [6]:
# Load participant information
partner_dict = {'sub-1001':'sub-2001','sub-2001':'sub-1001',
'sub-1002':'sub-2002','sub-2002':'sub-1002',
'sub-1003':'sub-2003','sub-2003':'sub-1003',
'sub-1004':'sub-2004','sub-2004':'sub-1004',
'sub-1005':'sub-2005','sub-2005':'sub-1005',
'sub-1006':'sub-2006','sub-2006':'sub-1006',
'sub-1007':'sub-2007','sub-2007':'sub-1007',
'sub-1008':'sub-2008','sub-2008':'sub-1008',
'sub-1009':'sub-2009','sub-2009':'sub-1009',
'sub-1010':'sub-2010','sub-2010':'sub-1010',
'sub-1011':'sub-2011','sub-2011':'sub-1011',
'sub-1012':'sub-2012','sub-2012':'sub-1012',
'sub-1013':'sub-2013','sub-2013':'sub-1013',
'sub-1015':'sub-2015','sub-2015':'sub-1015'}
dyad_ids = [1, 1, 2, 2, 3, 3, 4, 4,
5, 5, 6, 6, 7, 7, 8, 8,
9, 9, 10, 10, 11, 11,
12, 12, 13, 13, 14, 14]
subject_ids = ['sub-1001','sub-2001',
'sub-1002','sub-2002',
'sub-1003','sub-2003',
'sub-1004','sub-2004',
'sub-1005','sub-2005',
'sub-1006','sub-2006',
'sub-1007','sub-2007',
'sub-1008','sub-2008',
'sub-1009','sub-2009',
'sub-1010','sub-2010',
'sub-1011','sub-2011',
'sub-1012','sub-2012',
'sub-1013','sub-2013',
'sub-1015','sub-2015']
nsubj = len(subject_ids)
In [7]:
# create friend matrix and nonfriend_matrix
friend_matrix = np.zeros((len(subject_ids), len(subject_ids)))
friend_matrix[:] = np.nan
nonfriend_matrix = np.zeros((len(subject_ids), len(subject_ids)))
nonfriend_matrix[:] = np.nan
for i, sub_i in enumerate(subject_ids):
for j, sub_j in enumerate(subject_ids):
if (sub_i == partner_dict[sub_j]) or (partner_dict[sub_i] == sub_j):
friend_matrix[i, j] = 1
friend_matrix[j, i] = 1
else:
nonfriend_matrix[i, j] = 1
nonfriend_matrix[j, i] = 1
# also exclude self
if i == j:
nonfriend_matrix[i, j] = np.nan
friend_matrix[i, j] = np.nan
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))
sns.heatmap(friend_matrix, xticklabels=subject_ids, yticklabels=subject_ids, ax=axs[0], cbar=False, cmap='RdBu')
axs[0].set_title('Friend Pairs')
sns.heatmap(nonfriend_matrix, xticklabels=subject_ids, yticklabels=subject_ids, ax=axs[1], cbar=False, cmap='RdBu')
axs[1].set_title('Non-friend Pairs')
plt.show()
In [8]:
# Load dictionary of dictionary of searchlight objects contains RDMs based on the cross-validated Euclidean distance between the neural patterns across trials and runs
# `dyad_rdms` is formatted like: dyad_rdms[sub_id][cond_id]
# e.g., {'sub-1001' : {'Meta': searchlight_object, 'Other': searchlight_object},
# 'sub-1002' : {'Meta': searchlight_object, 'Other': searchlight_object},
# ...}
# pickled file is available at https://osf.io/t5cmd/
try:
with open('../data/dyad_rdms.pickle', 'rb') as handle:
dyad_rdms = pickle.load(handle)
except:
# if not downloaded, load from OSF
dyad_rdms = pickle.load(urlopen("https://osf.io/t5cmd/data/dyad_rdms.pickle", 'rb'))
In [9]:
# Load mask
mask_fn = 'MNI152NLin2009c-Schaefer2018-DMN21-k100_3mm.nii.gz'
mask_file = nib.load(mask_fn)
mask_img = mask_file.get_fdata()
mask = mask_img==1
x, y, z = mask_img.shape
In [ ]:
# Run DyRSA for every possible pair, then can filter friend vs non-friend pairs later
# first grab nvoxels from first subject
voxel_indices = list(dyad_rdms['sub-1001']['Meta'].rdm_descriptors['voxel_index'])
nvoxels = len(voxel_indices)
# initialize result matrices
pairwise_zmat_brain = np.empty((nvoxels, nsubj, nsubj))
pairwise_zmat_brain[:] = np.nan
# loop through all pairs and compute DyRSA
for subject_index, subject_id in tqdm(enumerate(subject_ids)):
for subject_index2, subject_id2 in enumerate(subject_ids):
if subject_index == subject_index2:
pairwise_zmat_brain[:, subject_index, subject_index2] = np.nan
else:
# Neural Meta(A(B(A))) ~ Neural Other(B(A))
subj_sl_RDM = dyad_rdms[subject_id]['Meta']
partner_sl_RDM = dyad_rdms[subject_id2]['Other'] #can also replace with other model RDMs
slresults = Parallel(n_jobs=-1)(delayed(spearmanr)(partner_sl_RDM.dissimilarities[x_index],
x.dissimilarities[0,]) for x_index, x in enumerate(subj_sl_RDM))
spearmanrhos_brain = [float(e[0]) for e in slresults]
pairwise_zmat_brain[:, subject_index, subject_index2] = np.array([np.arctanh(r) for r in deepcopy(spearmanrhos_brain)])
In [20]:
# remove friend values and flatten across voxels to shape nvox, sum(nonfriend_matrix)=728
all_nonfriend_pairs_flat = np.zeros((pairwise_zmat_brain.shape[0],
int(np.sum(nonfriend_matrix==1))))
# remove non-friend values and flatten across voxels to shape nvox, nsubj
all_friend_pairs_flat = np.zeros((pairwise_zmat_brain.shape[0], nsubj))
# loop through voxels
for v in range(pairwise_zmat_brain.shape[0]):
all_nonfriend_pairs_flat[v] = pairwise_zmat_brain[v,:,:][nonfriend_matrix==1]
all_friend_pairs_flat[v] = pairwise_zmat_brain[v,:,:][friend_matrix==1]
In [ ]:
# compute the mean of all_friend_pairs_flat
actual_group_mean = np.nanmean(all_friend_pairs_flat, axis=1) # mean at every voxel
# generate null distribution of bootstrap samples
nboot = 5000
zstat, perm_p = dyrsa.run_permutations(nsubj, nboot,
actual_group_mean,
all_nonfriend_pairs_flat,
replacement=True,
write_zstat=True)
In [ ]:
# Write unthresholded result
results_brain = np.zeros((x*y*z))
results_brain[voxel_indices] = zstat
results_brain = results_brain.reshape((x, y, z))
unthresh_img = new_img_like(mask_file, results_brain)
In [45]:
# plot
view_img_on_surf(unthresh_img,threshold=1)
Out[45]:
In [ ]:
# Apply FDR correction and write thresholded result
mask_pval = [pv for mask_ind, pv in zip(mask.ravel(), perm_p) if mask_ind]
fdr_thresh = fdr(np.array(mask_pval), q=.05)
results_brain2 = np.zeros((x*y*z))
results_brain2[voxel_indices] = np.array([m if p<fdr_thresh else np.nan for m, p in zip(zstat, perm_p)])
results_brain2 = results_brain2.reshape((x, y, z))
thresh_img = new_img_like(mask_file, results_brain2)
In [43]:
# plot
view_img_on_surf(thresh_img,threshold=.001)
Out[43]:
In [ ]: